# coding:utf-8
# Author : hiicy redldw
# Date : 2019/02/19
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
# plt.rcParams['figure.figsize'] = (20,10)
# 随机行走时间序列
def random_walk():
    np.random.seed(1291)
    z=np.random.normal(0.1,2,100) # 均值，标准差 100个数
    y=np.cumsum(z)
    # print(z.std())
    # print(z)
    fig,ax1 = plt.subplots()
    plt.plot(z,label='white noise')
    plt.plot(y,label='Random Walk')
    plt.legend()
    plt.show()
    mean1=np.round(np.mean(y[:20]),4);mean2=np.round(np.mean(y[-20:]),4)
    std1=np.round(np.std(y[:20]),4);std2=np.round(np.std(y[-20:]),4)

    print(f'前20个数据点的均值为{mean1:4f},标准差为{std1:4f}')
    print('\\')
    print(f'后20个数据点的均值为{mean2:4f},标准差为{std2:4f}')

# random_walk()

def changjiang():
    path = r'F:\Resources\Dataset\monthly-flows-chang-jiang-at-han.csv'
    parser = lambda date:pd.datetime.strptime(date,"%Y-%m")
    df1 = pd.read_csv(path,engine='python',skipfooter=3,names=['YearMonth','WaterFlow'],
                      parse_dates=[0], infer_datetime_format=True, date_parser=parser, header=0)
    df1.YearMonth = pd.to_datetime(df1.YearMonth)
    print(df1.YearMonth)
    df1.set_index("YearMonth",inplace=True) # 把标签设为索引
    df1.plot()
    # 该数据具备很强的不同长度的周期性
    plt.show()
# changjiang()


def airline_passengers():#该数据包含极强的趋势要素和周期要素
    path = r'F:\Resources\Dataset\international-airline-passengers.csv'
    parser = lambda date:pd.datetime.strptime(date,"%y-%b")
    df2 = pd.read_csv(path,engine='python',skipfooter=3,names=['YearMonth','Passenger'],
                      header=0)
    # print(df2.YearMonth)
    # df2.YearMonth = df2.YearMonth.str[:4]+'19'+df2.YearMonth.str[-2:]
    df2.YearMonth = pd.to_datetime(df2.YearMonth,infer_datetime_format=True)
    df2.set_index("YearMonth",inplace=True)
    # 先确定差分阶数
    order = 1
    diff1 = df2.Passenger.diff(order)[1:]
    # print(len(df2),diff1,len(diff1))
    logdiff1 = np.log(df2.Passenger).diff(order)[order:]
    # print(len(diff1),logdiff1,len(logdiff1))
    adftest = sm.tsa.stattools.adfuller(diff1)
    adftestlog = sm.tsa.stattools.adfuller(logdiff1)
    # 以p=0.05 作为界值
    print('ADF test result on Difference shows test statistic is {0[0]:f} '
          'and p-value is {0[1]:f}'.format(adftest[:2]))
    print('ADF test result on Log Difference shows test statistic is {0[0]:f}'
          'and p-value is {0[1]:f}'.format(adftestlog[:2]))
    import warnings
    with warnings.catch_warnings():
        warnings.filterwarnings('ignore')
        kpsstest = sm.tsa.stattools.kpss(diff1)
        kpsstestlog = sm.tsa.stattools.kpss(logdiff1)
    print('KPSS test result on Difference shows test statistic is %f \
            adn p-value is %f' % (kpsstest[:2]))
    print('KPSS test result on Log Difference shows test statistic is {0[0]:f}'
          'and p-value is {0[1]:f}'.format(kpsstestlog[:2]))
    # ACF和PACF情况
    print(diff1)
    fig,ax = plt.subplots()
    ax1 = fig.add_subplot(221)
    sm.graphics.tsa.plot_acf(diff1,ax=ax1,lags=50)
    ax2 = fig.add_subplot(222)
    sm.graphics.tsa.plot_pacf(diff1,ax=ax2)

    ax3 = fig.add_subplot(223)
    sm.graphics.tsa.plot_acf(logdiff1, ax=ax3)
    ax4 = fig.add_subplot(224)

    sm.graphics.tsa.plot_pacf(logdiff1, ax=ax4)

    # print(df2.head())
    # logdiff1.plot()
    plt.show()

# airline_passengers()
import keras.models as KModels
import keras.layers as KLayers
from scipy.signal import periodogram
import warnings
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import mean_squared_error


def test_stationarity(timeseries,window=12):
    """
    测试平稳性
    """
    df = pd.DataFrame(timeseries)
    # 观测移动平均和移动均方差随时间的变化图
    df['Rolling.Mean'] = timeseries.rolling(window=window).mean()
    df['Rolling.Std'] = timeseries.rolling(window=window).std()
    adftest = sm.tsa.stattools.adfuller(timeseries)
    adfoutput=pd.Series(adftest[0:4],index=['统计量','p-值 ','滞后量','观测值数量'])
    for key,value in adftest[4].items():
        adfoutput['临界值 {%s}' % key] = value
    return (adfoutput,df)

def CalculateCycle(ts,lags=36):
    import peakutils as peak
    from statsmodels.tsa.stattools import acf
    from scipy import signal
    acf_x,acf_ci = acf(ts,alpha=0.05,nlags=lags)
    fs=1
    # 周期图法 识别周期
    f,Pxx_den = signal.periodogram(acf_x,fs)
    index=peak.indexes(Pxx_den)
    # print('index',index)
    cycle=(1/f[index[0]]).astype(int)
    fig = plt.figure()
    ax0 = fig.add_subplot(111)
    plt.vlines(f,0,Pxx_den)
    plt.plot(f,Pxx_den,marker='o',linestyle='none',color='red')
    plt.title('Identified Cycle of %i' % (cycle))
    plt.xlabel('frequency [Hz]')
    plt.ylabel('PSD [V**2/Hz]')
    plt.show()
    return (index,f,Pxx_den)


def changjian_arima():
    path = r'F:\Resources\Dataset\monthly-flows-chang-jiang-at-han.csv'
    parser = lambda date: pd.datetime.strptime(date, "%Y-%m")
    df1 = pd.read_csv(path, engine='python', skipfooter=3, names=['YearMonth', 'WaterFlow'],
                      parse_dates=[0], infer_datetime_format=True, date_parser=parser, header=0)

    df1.YearMonth = pd.to_datetime(df1.YearMonth)
    # print(df1)
    df1.set_index("YearMonth", inplace=True)  # 把标签设为索引
    # m = CalculateCycle(df1)

    cutoff = 24
    train:pd.Series = df1.WaterFlow[:-cutoff]
    test = df1.WaterFlow[-cutoff:]
    Seasonality=12

    # 以12差分 去掉季节性
    waterFlowS12=train.diff(Seasonality)[Seasonality:]
    adftestS12=sm.tsa.stattools.adfuller(waterFlowS12)
    print(adftestS12)
    print('ADF test result on Difference shows test statistic is {0[0]:f} '
          'and p-value is {0[1]:f}'.format(adftestS12[:2]))
    nlag=36
    xvalues=np.arange(nlag+1)
    acfS12,confiS12=sm.tsa.stattools.acf(waterFlowS12,nlags=nlag,alpha=0.05,fft=False)
    confiS12 = confiS12 + confiS12.mean(1)[:,None]  # 往里加方号，
    print('confiS12:{conf},   {conf.shape}'.format(conf=confiS12))

    # fig = plt.figure()
    # ax0 = fig.add_subplot(221)
    # waterFlowS12.plot(ax=ax0)

    # ax1=fig.add_subplot(222)
    #  12的倍数滞后项都没有显著的数据值 所以12的周期是可以的
    # sm.graphics.tsa.plot_acf(waterFlowS12,lags=nlag,ax=ax1)

    fig = plt.figure()
    waterFlowS12d1=train.diff(1)[1:]
    print('waterflows12d1:',waterFlowS12d1)
    ax0 = fig.add_subplot(221)
    # sm.graphics.tsa.plot_acf(waterFlowS12d1,ax=ax0,lags=48)

    ax1 = fig.add_subplot(222)
    # sm.graphics.tsa.plot_pacf(waterFlowS12d1,ax=ax1,lags=48)
    # ACF图表明可能需要1个季节性预测误差滞后项，即差分滞后项
    # 一个季节性ARIMA模型order(p,d,q)
    mod1 = sm.tsa.statespace.SARIMAX(train,order=(0,1,0),
                                     seasonal_order=(0,1,1,12)).fit()
    mod2 = sm.tsa.statespace.SARIMAX(train, order=(1, 1, 1),
                                     seasonal_order=(0, 1, 1, 12)).fit()
    pred=mod1.predict()
    print(mod1.summary())
    subtrain = train['1960':'1970']
    MAPE = (np.abs(train - pred)/train).mean()
    subMAPE = (np.abs(subtrain-pred['1960':'1970'])/train).mean()
    fig = plt.figure()
    ax0 = fig.add_subplot(211)
    plt.plot(pred,label='Fitted')
    plt.plot(train,color='red',label='Original')
    plt.legend(loc='best')
    plt.title('SARIMA(0,1,0) (0,1,1,12) Model,MAPE={:.4f}'.format(MAPE))

    ax1 = fig.add_subplot(212)
    plt.plot(pred['1960':'1970'], label='Fitted')
    plt.plot(subtrain, color='red', label='Original')
    plt.legend(loc='best')
    plt.title('Details from 1960 to 1970, MAPE = {:.4f}'.format(subMAPE))
    # 进行预测
    fig1 = plt.figure()
    forecast1 = mod1.predict(start='1976-12-01',end='1978',dynamic=True)
    forecast2 = mod2.predict(start='1976-12-01',end='1978',dynamic=True)
    MAPE1 = ((test-forecast1).abs()/test).mean()*100
    MAPE2 = ((test-forecast2).abs()/test).mean()*100
    plt.plot(test,color='black',label='Original')
    plt.plot(forecast1,color='green',label='Model 1 : SARIMA(0,1,0)(0,1,1,12)')
    plt.plot(forecast2,color='red',label='Model 2 : SARIMA(1,1,1)(0,1,1,12)')
    plt.legend(loc='best')
    plt.title('Model 1 MAPE=%.f%%; Model 2 MAPE=%.f%%' % (MAPE1,MAPE2))



    # print(train,type(train))
    # fig = plt.figure()
    # ax0 = fig.add_subplot(221)
    # adftest, dftest0 = test_stationarity(train)
    # dftest0.plot(ax=ax0)
    # # print('原始数据平稳性校验')
    # # print(adftest)


    # ax1=fig.add_subplot(222)
    # adftest,dftest1 = test_stationarity(train['1960':'1975'])
    # dftest1.plot(ax=ax1)
    # print('局部数据平稳性校验')
    # print(adftest)
    plt.show()

    # 移动方差波动幅度大，非常强的季节性
def changjiang_lstm():
    def create_dataset(dataset,timestep=1,look_back=1,look_ahead=1):
        from statsmodels.tsa.tsatools import lagmat
        # print(dataset)
        ds = dataset.reshape(-1,1) # 给定一列，自动计算多少行
        dataX = lagmat(dataset,maxlag=look_back,trim='both',original='ex')
        # print(dataX)
        dataY = lagmat(dataset[look_back:],maxlag=look_ahead,trim='backward',original='ex')
        # print(dataX.shape,dataY.shape)
        # reshape and remove redundent rows
        print(dataX.shape)
        dataX = dataX.reshape(dataX.shape[0],timestep,dataX.shape[1])[:-(look_ahead-1)]
        print('dataX',dataX.shape,dataY[:(look_ahead+1)].shape)
        return np.array(dataX),np.array(dataY[:-(look_ahead-1)])

    path = r'F:\Resources\Dataset\monthly-flows-chang-jiang-at-han.csv'
    parser = lambda date: pd.datetime.strptime(date, "%Y-%m")
    df1 = pd.read_csv(path, engine='python', skipfooter=3, names=['YearMonth', 'WaterFlow'],
                      parse_dates=[0], infer_datetime_format=True, date_parser=parser, header=0)

    df1.YearMonth = pd.to_datetime(df1.YearMonth)
    # print(df1)
    df1.set_index("YearMonth", inplace=True)  # 把标签设为索引
    # m = CalculateCycle(df1)
    from sklearn import preprocessing

    cutoff = 24
    train: pd.Series = df1.WaterFlow[:-cutoff]
    test = df1.WaterFlow[-cutoff:]
    scaler = preprocessing.MinMaxScaler(feature_range=(0,1)) #标准化 【0,1】
    trainstd = scaler.fit_transform(train.values.astype(float).reshape(-1,1))
    teststd = scaler.fit_transform(test.values.astype(float).reshape(-1,1))

    lookback=60
    lookahead=24
    timestep=1
    trainX,trainY = create_dataset(trainstd,timestep=timestep,look_back=lookback,look_ahead=lookahead)
    print('trainX',trainX.shape,'trainY',trainY.shape)
    batch_size = 1
    model = KModels.Sequential()

    model.add(KLayers.LSTM(48,batch_size=batch_size,input_shape=(1,lookback),kernel_initializer='he_uniform'))
    model.add(KLayers.Dense(lookahead))
    model.compile(loss='mean_squared_error',optimizer='adam')
    # model.fit(trainX,trainY,batch_size=1,epochs=20,verbose=0)
    print(df1.WaterFlow['1972':'1976'])
    feedData=scaler.transform(df1.WaterFlow['1972':'1976'].reshape(-1,1))
    feedX = (feedData).reshape(1,1,lookback)
    feedX = (feedX)
    prediction1 = model.predict(feedX)
    predictionRaw = scaler.inverse_transform(prediction1.reshape(-1,1))
    actual1=df1.WaterFlow['1977':'1978'].copy().reshape(-1, 1)
    MAPE = (np.abs(predictionRaw - actual1)/actual1).mean()
    plt.plot(predictionRaw,label='Prediction')
    plt.plot(actual1,label='Actual')
    plt.title('MAPE = %.4f' % MAPE)
    plt.legend(loc='best')
    plt.xlim((0,23))
    plt.xlabel('Month')
    plt.show()


"""
检验平稳性
"""
changjiang_lstm()




















